load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"


location_f="/Users/ivana/WORK_BSC/data_to_publish/CCSM4som/"

filsw= systemfunc ("ls "+location_f+"lowice*2dvars*.nc")
filsc= systemfunc ("ls "+location_f+"ctrl*2dvars*.nc")

fw    = addfiles (filsw, "r")
fc    = addfiles (filsc, "r")

ListSetType (fw, "join")
ListSetType (fc, "join")
print(fw)
print(fc)

;*******************************************;


icefracc = fc[:]->ICEFRAC  
icefracw = fw[:]->ICEFRAC

printVarSummary(icefracc)

icec =new((/6,12,96,144/),typeof(icefracc))
icew =new((/12,12,96,144/),typeof(icefracc))

do i=0,5
icec(i,:,:,:)=clmMonTLL(icefracc(i,:,:,:))
end do

do i=0,11
icew(i,:,:,:)=clmMonTLL(icefracw(i,:,:,:))
end do


lat=fc[0]->lat
lon=fc[0]->lon

nlat=dimsizes(lat)
nlon=dimsizes(lon)

re   = 6.37122e06

   rad  = 4.0 * atan(1.0) / 180.
   con  = re * rad                 
   clat = cos(lat * rad)           ; cosine of latitude
   
   dlon = (lon(2) - lon(1))        ; assume dlon is constant
   dlat = (lat(2) - lat(1))        ; assume dlat is constant

   dx   = con * dlon * clat        ; dx at each latitude
   dy   = con * dlat               ; dy is constant
   dydx=dx
   dydx = dy * dx                  ; dydx(nlat)
   
   wgt  = new((/nlat, nlon/), typeof(lat))
   wgt  = conform(wgt, dydx, 0)

icesumc = wgt_areasum2(icec(:,:,60:95,:),wgt(60:95,:),0)
icesumc = icesumc*10^(-12)

icesumw = wgt_areasum2(icew(:,:,60:95,:),wgt(60:95,:),0)
icesumw = icesumw*10^(-12)

month = ispan(0,11,1)*1.0
;****************************************
; Create plot
;****************************************
wks = gsn_open_wks("pdf","icearea_monthly_ccsm4som")

   plot    = new (1,"graphic")

   res                      = True          ; individual plot
   res@gsnDraw              = False
   res@gsnFrame             = False
   res@xyLineThicknessF     = 2.0
   res@xyMonoLineColor         =True
   res@xyLineColor         = "black"

   res@tmXBMode        = "Explicit"
   res@tmXBValues      = month

   res@tmXBLabels      = (/"J","F","M","A","M","J"\
              ,"J","A","S","O","N","D"/)

   res@xyDashPattern        = 0                  ; Make curves all solid

   res@vpHeightF= 0.4                    ; change aspect ratio of plot
   res@vpWidthF = 0.8

   res@trXMinF  =  0
   res@trXMaxF  = 11
   res@trYMinF=0
   res@trYMaxF=16

   res@tiXAxisFontHeightF = 0.010
   res@tiYAxisFontHeightF = 0.015

   res@gsnLeftString = ""
   res@tiXAxisString = ""
   res@tiYAxisString = ""

   icesumc!0 = "ncl_join"
   icesumc!1= "time"
   icesumc&time= month

   icesumw!0 = "ncl_join"
   icesumw!1= "time"
   icesumc&time= month


   plot(0)  = gsn_csm_xy (wks,month,icesumc,res)
   res@xyLineColor         = "cyan3"
   plot1  = gsn_csm_xy (wks,month,icesumw,res)

   overlay(plot(0), plot1)

   ;***** make panel *****


   resP                     = True          ; panel resources
   resP@gsnMaximize         = True

   resP@txFontHeightF = 0.015
   gsn_panel(wks,plot,(/1,1/),resP)
